Sourcing Data
The data for this analysis is sourced from the website of the Central Statistics Office - Ireland’s agency tasked with keeping track of all the numbers. A census was performed in 2016, and aside from the standard information based on where people reside, detailed information has been released that outlines where people work or attend school. Further information can be found at CSO.
library(dplyr)
Attaching package: 㤼㸱dplyr㤼㸲
The following objects are masked from 㤼㸱package:stats㤼㸲:
filter, lag
The following objects are masked from 㤼㸱package:base㤼㸲:
intersect, setdiff, setequal, union
library(magrittr)
library(sf)
Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3
library(downloader)
library(stringr)
library(purrr)
Attaching package: 㤼㸱purrr㤼㸲
The following object is masked from 㤼㸱package:magrittr㤼㸲:
set_names
if (!dir.exists("data")) {
dir.create("data")
}
if (!dir.exists("output")) {
dir.create("output")
}
if (!dir.exists("plots")) {
dir.create("plots")
}
if (!file.exists("data/Workplace_Zones_ITM.shp")) {
# download
downloader::download("http://www.cso.ie/censusfiles/Workplace_Zones_ITM.zip",
dest = "data/wz.zip", mode = "wb")
unzip("data/wz.zip", exdir = "data")
file.remove("data/wz.zip")
}
if (!file.exists("data/wz.xlsx")) {
# download
downloader::download("http://www.cso.ie/en/media/csoie/census/census2016/census2016results/saps/Workplace_zones_-_SAPS_2016.xlsx", dest = "data/wz.xlsx", mode = "wb")
}
if (!file.exists("data/wz_lookup.xlsx")) {
# download
downloader::download("http://www.cso.ie/en/media/csoie/census/census2016/census2016results/saps/Workplace_Zones_SAPs_Theme_breakdown.xlsx", dest = "data/wz_lookup.xlsx", mode = "wb")
}
wp_zones <- st_read("data/Workplace_Zones_ITM.shp")
Reading layer `Workplace_Zones_ITM' from data source `C:\Users\Admin\Documents\R\cso_bike\data\Workplace_Zones_ITM.shp' using driver `ESRI Shapefile'
Simple feature collection with 7219 features and 8 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 417471.5 ymin: 519663.7 xmax: 734481.1 ymax: 966896.3
epsg (SRID): NA
proj4string: +proj=tmerc +lat_0=53.5 +lon_0=-8 +k=0.99982 +x_0=600000 +y_0=750000 +ellps=GRS80 +units=m +no_defs
wp_data <- readxl::read_xlsx("data/wz.xlsx")
wp_lookup <- readxl::read_xlsx("data/wz_lookup.xlsx")
wp_lookup <- wp_lookup[-1,-c(1:2,5)]
The entire data for the Republic of Ireland is now downloaded. We need to subset the data for the area we are interested in. This broadly relates to the area within Dublin City Council. The statistics included in this data are reasonably vague - part of this is because the Census is required to be anonymous, but there is also ambigous aspects to how the people working, at school or at home are counted. Again, read the notes in the CSO description of this data to get more insight. I have done my best to intrepret it, and where possible I’ve explained my reasoning..
wp_zones %>% filter(COUNTY == 'DC') -> dub_wz
st_crs(dub_wz) <- 2157
# now lets join the polygons with the data
dub_wz <- dub_wz %>% left_join(wp_data)
Joining, by = "GUID"
Column `GUID` joining factor and character vector, coercing into character vector
# remove spurious data
dub_wz <- dub_wz[,-c(1:7)]
# modal share for transport
# subset fields
dub_wz %>% select(T11_C1, starts_with('T2')) -> dub_mode
dub_mode[is.na(dub_mode)] <- 0
dub_mode %<>% mutate(prop_cycling = T2_M2 / T2_T)
dub_mode %<>% mutate(prop_walking = T2_M1 / T2_T)
dub_mode %<>% mutate(prop_bus = T2_M3 / T2_T)
dub_mode %<>% mutate(prop_train = T2_M4 / T2_T)
dub_mode %<>% mutate(prop_motorcycle = T2_M5 / T2_T)
dub_mode %<>% mutate(prop_car = (T2_M5 + T2_M6 + T2_M7) / T2_T)
dub_mode %<>% mutate(prop_other_ns = (T2_M9 + T2_NS) / T2_T)
dub_mode_prop <- dub_mode %>% st_set_geometry(NULL) %>% round(2)
dub_mode_prop$guid <-
# calculate prop for journey time
dub_wz %>% select(starts_with('T4')) -> dub_journey_time
Warning in cbind(x[0:(framecol - 1)], cols) :
number of rows of result is not a multiple of vector length (arg 2)
# calculate time leaving home
dub_wz %>% select(starts_with('T5')) -> dub_leave_home
We now want to have a look at our data. At this moment in time, the best way to view R geospatial data is via the Leaflet web-mapping package, which is a Java based infrastructure. It allows interactive and close examination of the data in a way that would be familiar to anyone using maps on the internet. The geom_sf functions of the ggplot2 universe have not yet reached production standard, and also I had trouble installing them on my machine for some reason….
#first transform our wz data to wgs84, a web friendly projection that'll play nicely with leaflet
dub_mode_wgs84 <- dub_mode %>% st_transform(4326)
library(leaflet)
package 㤼㸱leaflet㤼㸲 was built under R version 3.4.4
#show a map of dublin city centre
d <- leaflet() %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$OpenMapSurfer.Grayscale, group = "Gray") %>%
addProviderTiles(providers$Thunderforest.OpenCycleMap, group = "CycleMap") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite")
d %>% addMarkers(lng= -6.271394, lat=53.345522, popup="The Dublin Quays")
#now we'd like to show the proportional breakdown of cycling mode
# first set up a colour palette - see https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf
# this also very useful - https://rstudio.github.io/leaflet/choropleths.html
pal_cycling <- colorBin(
palette = "YlOrRd",
domain = dub_mode_wgs84$prop_cycling)
d_wz_cycling <- d %>% addPolygons(data = dub_mode_wgs84,
weight = 0,
opacity = 1,
color = 'navy',
fillOpacity = 0.6,
fillColor = ~pal_cycling(prop_cycling))
d_wz_cycling %>% addLegend(pal = pal_cycling, values = dub_mode_wgs84$prop_cycling, position = "bottomright")
pal_car <- colorNumeric(
palette = "Reds",
domain = dub_mode_wgs84$prop_car)
d_wz_car <- d %>% addPolygons(data = dub_mode_wgs84,
weight = 0,
opacity = 1,
color = 'red',
fillOpacity = 0.6,
fillColor = ~pal_car(prop_car))
d_wz_car %>% addLegend(pal = pal_car, values = dub_mode_wgs84$prop_car, position = "bottomright")
# create a colour scheme for proportionality
bins <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
pal_prop <- colorBin("YlOrRd", domain = c(0, 1), bins = bins)
# add labels to each polygon
car_labels <- sprintf(
"%g Motorcycle<br/>%g Car Driver<br/>%g Car Passenger<br/><strong>%g Total</strong>",
dub_mode_wgs84$T2_M5, dub_mode_wgs84$T2_M6, dub_mode_wgs84$T2_M7, dub_mode_wgs84$T2_T
) %>% lapply(htmltools::HTML)
bike_labels <- sprintf(
"%g Bicycle<br/><strong>%g Total</strong>",
dub_mode_wgs84$T2_M2, dub_mode_wgs84$T2_T
) %>% lapply(htmltools::HTML)
d_wz_car_cycling <- d %>%
addPolygons(data = dub_mode_wgs84,
group = "Bicycle",
weight = 0,
opacity = 1,
color = 'navy',
fillOpacity = 0.6,
fillColor = ~pal_prop(prop_cycling),
label = bike_labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")
) %>%
addPolygons(data = dub_mode_wgs84,
group = "Car (Driver and Passenger)",
weight = 0,
opacity = 1,
color = 'red',
fillOpacity = 0.6,
fillColor = ~pal_prop(prop_car),
label = car_labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")
)
d_wz_car_cycling %<>% addLayersControl(
baseGroups = c("OSM (default)", "Gray", "CycleMap", "Satellite"),
overlayGroups = c("Bicycle", "Car (Driver and Passenger)"),
options = layersControlOptions(collapsed = FALSE)
) %>%
addLegend(pal = pal_prop, values = bins, position = "bottomright", title = "Proportion of Total Daytime Population")
d_wz_car_cycling
NExt on the list: make this map beautiful..! make it useful - able to view cycling proportion with comparative shares of other modes - see what effect a good transport link has on the outcome - see what the distribution of travel durations is